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Unci as 


INTRODUCTION 


The applicants have developed a finite element based approach for the solution of 3D 
compressible flows. The procedure enables flow solutions to be obtained on 
letrahedral discretisations of computational domains of complex form. A further 
development has been the incorporation of a solution adaptive mesh strategy in which 
the adaptivity is achieved by complete remeshing of the solution domain. During the 
previous year the applicants have been working with the Advanced Aerodynamics 
Concepts Branch at NASA Ames Research Center with an implementation of the basic 
meshing and solution procedure. The objective of the work to be performed over this 
twelve month period was the transfer of the adaptive mesh technology and also the 
undertaking of basic research into alternative flow algorithms for the Euler 
equations on unstructured meshes. 

COMPUTER SOFTWARE 

J. Peiro visited Ames Research Center during the period 30th April - 4 May, 1990. 
He undertook the delivery of new versions of the Taylor-Galerkin flow solver 
(ARCTG), the surface triangulator (ARCST) and the volume generator (ARCVT). In 
addition he answered queries about the basic algorithms and code organisation and 
provided information on the modifications present in the updated versions of the code. 
As a major task under the current grant, the computer code necessary to perform the 
adaptive remeshing was provided. This code (ARCERR) provides new mesh parameter 
information based upon an error estimation procedure which uses the second 
derivatives of the solution. Dr. Erickson of NASA Ames R.C. requested additional codes 
and we agreed to make these available. These codes provided the ability to generate 3D 
meshes from a 2D unstructured triangular grid (ARCGB) and for graphic post- 
processing of the geometry and flow solution (ARCUT). To assist the group at NASA 
Ames R.C. with the process of familiarisation with the adaptive remeshing codes, J. 
Peiro undertook in London the adaptive solution of the problem of flow at Mach 1 .75 


2 



past a 6 degree double wedge aerofoil. This work, which is attached here in the form 
of Appendix A, was written up in the form of a report and copies provided to Dr. 
Erickson. 


BASIC RESEARCH ON FLOW ALGORITHM DEVELOPMENT 


Theory 

As a continuation of the work undertaken last year, further investigative work has 
been undertaken into the possibility of developing a least-squares based Euler solver 
for implementation on unstructured grids. The current generation of solvers are 
considered to be far too sensitive to grid quality. The driving force behind these 
investigations continues to be the hope that a successful implementation will allow 
for the incorporation of quadratic elements, which should prove to be less sensitive 
to the quality of the grid being employed. 


The formulation which has been tested works with the governing equation 


auaFdG au A au Q au 

9t + 3x + 3y “ 3t + 9x + 3y 


0 


( 1 ) 


where A and B are the standard jacobian matrices and U is the vector of the 
conservation variables. This equation is discretised in time by the using a backward 
difference approximation, which can be written as 

AU n 3(AU) n 9(AU) [3F dG_T (2) 

At +A 3x +B ay = 'Lax + 3yJ [d) 

where AU = U n+1 - U n is the increment in U over the timestep At and A n = A(U n ), 
B n * B(U n ). The least squares method is based upon the minimisation of the L 2 norm 
of this equation with respect to U n+1 . If the domain is divided into n-noded triangular 
elements (with n=3 for linear elements and n=6 for quadratic elements), then 
approximations for All, F n and G n can be taken in the form 
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AU = AUjNj 


(3) 


F n = F, n Ni G n = G, n Ni 

where N ( denotes the finite element shape function associated with node i of the mesh. 
A piecewise constant representation is employed for the matrices A n and B n . Upon 
insertion of these expansions into equation ( 2 ), the least squares requirement leads 
to 

Kij AUj = fj (4) 


where the entries in the matrix K and the right hand side vector f are defined by 


Kij = J [INj + At A n ^ + AtB n ^ i ] T [IN j +AtA n ^ + AtB n f^]dQ 


fi 


3F n 3G n 


.JnN 1 **.*"f + 4 ,B"^ 1 T | ^ r+ « rl dD 

When steady state conditions are reached, equation (4) reduces to 


(5) 


J 


3F n 3G n 


< N 'l *1 


. , r » n 3Nj Al on 3Nj T n 3U n B n 3U n . . HO . 
+ At [A 3 ^- + At B ] 1 [A + At B a „ ] } dQ - 0 


3y 


3y 


( 6 ) 


and the first term here can be seen to represent a Galerkin form of the steady state 
Euler equation set, while the second term represents a added diffusion. It can be 
observed that the amount of added diffusion depends upon the magnitude of the time 
step and on the jacobian matrices. 
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The time step employed can be taken to be uniform everywhere or, for the steady 
computations which are of interest here, leal timesteps can be employed. It has been 
found that the convergence rate to steady state is enhanced by the use of local 
timestepping but that it also leads to the creation of oscillations in the solution in the 
vicinity of high gradient regions, especially when adaptive mesh refinement is 
employed. It should also be noted that numerical integration is required for the 
evaluation of the integrals appearing in equation (5) when curved-sided triangular 
elements are used in the discretisation of the domain. The solution to equation (4) is 
achieved by the incomplete Choleski conjugate gradient method. 

An Adaptive Refinement Strategy 

The diffusion term in equation (6) can be identified with the L 2 norm of the residual 
and so provides a natural method for indicating where the mesh resolution is 
insufficient. The standard adaptive remeshing procedure has been modified to use this 
term as an error indicator. 

Numerical Examples 

Several examples, involving flows with a supersonic freestream, have been solved 
with the above algorithm using both linear and quadratic elements. It can be 
generally seen from the results presented here that quadratic meshes give more 
accurate results, but they do suffer from slower convergence and tend to exhibit 
some undesirable oscillations. The flow tangency condition which should be satisfied 
at solid walls is imposed weakly via a least squares functional term. A characteristic 
analysis is carried out along the computational boundaries. 

(a) Double Wedge Aerofoil. Details of the meshes employed in the solution of this 
example are shown in figure 1 . The free stream Mach number was 2. Also shown in 
figure 1 are the convergence rates achieved and the distribution of the pressure over 
the wedge. 
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(b) 20° Ramp. The free stream Mach number for these computations is 3 and the 
results are shown in figure 2. The application of the adaptive mesh strategy to this 
problem is illustrated in figure 3. 

(c) Double Compression Corner. The free stream Mach number for this example is 
2 and the results obtained are shown in figure 4. The superiority of the quadratic 
elements over linear elements is most apparent in this example. The solution 
obtained by adaptive remeshing is shown in figure 5. 

(d) Supersonic Flow Past a Cylinder. Supersonic flow at a free stream Mach number 
of 3 past a circular cylinder is considered in this example. Only one quarter of the 
cylinder is considered as shown in figure 6. Since the implicit artificial diffusion of 
the scheme is directly related to the jacobian matrices, it is apparent that problems 
may be expected in the vicinity of stagnation points. In this case, a small amount of 
explicit artificial viscosity had to be added to maintain the stability of the scheme. An 
adapted mesh solution is shown in figure 7. 

Current Work 

The method developed to date has been shown to be capable of producing encouraging 
results for compressible flows at low supersonic Mach numbers. However, 
numerical experiments with quadratic elements have indicated a tendency to produce 
oscillatory solutions. In addition, when refinement is employed, the method tends to 
generate oscillations in the vicinity of regions of steep gradient. These problems can 
be overcome by using a small amount of explicitly added artificial diffusion and this 
is relatively straightforward in the case of linear elements as several such methods 
already exist. Extending these methods to produce a diffusion operator on quadratic 
elements which is consistent with the algorithm is a task which is still to be 
performed. 
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CONCLUSIONS 


Updated versions of a finite element based set of codes for the solution of the 
compressible Euler equations on general unstructured meshes have been successfully 
implemented at NASA Ames R.C. In addition, a facility for obtaining mesh adaptation 
by complete remeshing has been provided. Work on algorithmic development, with 
the objective of improving the solution quality on general grids, has proved of 
limited success and it is apparent that further efforts are needed to improve the 
situation in this area. 


PUBLICATION 

J. Peiro, J. Peraire and K. Morgan, 'The computation of 3D aerodynamical flows 
using unstructured meshes’, in Science and Engineering on Supercomputers, Edited 
by Eric J. Pitcher, CML Publications and Springer Verlag, 103-111, 1990. 


7 












Pressure contours 



First adaptation 1658 linear elements 




Second adaptation 1092 linear elements 



First adaptation 523 quadratic elements 


FIGURE 3 Mach 3 flow over a 20° ramp - adapted mesh solutions 
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FIGURE 5 Mach 2 flow over a double compression corner - adapted mesh solution 
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Pressure contours 


FIGURE 7 Mach 3 flow over a cylinder - adapted mesh solutions 
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Report on the solution of the 
mach 1.75 double wedge 


This report presents the results obtained for the solution of the flow 
past a 6° double wedge at M«=1.75 and a=0°. These results have been 
produced with the same codes currently in use in the Ames Research 
Center. The source files can be found in the SWANSEA account in 
RALPH and the names used are: 


• arcst : 

• arcvt : 

• setboun 

• arctg : 

• arcgb : 

• arcerr : 

• arccut : 


Surface triangulator 

Tetrahedral mesh generator 

Sets the boundary conditions for the 3D mesh 

3D Taylor-Galerkin Euler solver 

Generator of initial 3D background meshes from 

2D meshes 

Program for error analysis and definition of 
background meshes from computed solutions 
Plotting program 


The definition of the boundary of the computational domain 
employed in our computations was provided by Jahed. The data file 
(dwlp.dat) is enclosed in pages 4-5. A sketch of the computational 
domain is also included in page 6. 

The background grid for the generation of the first mesh has been 
obtained from a 2D mesh using the program arcgb. The input data 
file (dwlp.bck) and a drawing of the 2D mesh are displayed in pages 
7-8. The final background grid data (dwlp.bac) can be found in pages 
9 - 12 . 


The generated mesh, which will referred to as first mesh, contains 
42905 elements, 10180 nodes and 11230 boundary faces. Page 13 shows 
a view of the surface mesh in the plane y=0.1. A set of preliminary 
runs was performed on this mesh in order to assess the best choice of 
parameters for the control file to be used in the code arctg. We 
selected a diffusion coefficient of 1.1 and a smoothing type 2 and tried 
different types of boundary conditions (parameter istron): weak (0), 
strong (2) and mixed (1). The computed distributions of the pressure 
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coefficient on the surface of the wedge at y=0.05 are displayed in 
page 14. 

The solution obtained presents oscillations in both shocks and also in 
the expansion fan. The use of strong boundary conditions produces 
the best results giving the smallest oscilations in the trailing edge 
and the expansion fan. Since all the options produced oscillations on 
the leading edge shock we opted for the use of a bigger pressure 
switch coefficient. The control file (dwlp.con) employed for the final 
run is displayed in page 15. The flow solution obtained after 4000 
timesteps starting from free stream values is displayed as follows: 

page 16: pressure contours on the plane y=0.1 

page 17: mach number contours on the plane y=0.1 

page 18: pressure coefficient on the su rface of the 

wedge at y=0.05 compared with the exact 
solution 

page 19: pressure coefficient at the line z=0.0 y=0.1 

Note that the wiggles at the trailing edge and the expansion fan have 
been substantially reduced. The leading edge oscillations, although 
smaller, still are present. The magnitude of the oscillations will, 
however, decrease after the remeshing. 

This solution has been utilised to produce a second adapted mesh. The 
code used is arcerr and the "key" variable for the error indication is 
the mach number. The values employed to produce the background 
mesh for the generation of the second mesh are: 

^scaling = 0.0015 
^minimum = 0.009 
^maximum = 0.2 
stretching = 1.0 

The criterion for chosing these values was: 

1. - Use a smaller 6 minimum than in the first mesh 

2. - Keep the same value 6 maximum 

3. - Choose 5 sca ]ing such that the number of elements to be 
generated is aproximately equal or slightly higher than the one used 
in the first mesh (42905 elements). 

The estimated number of elements given by arcerr was 52450. The 
actual number of elements generated by arcst and arcvt in the 
second adapted mesh is 69742. The generated surface mesh in the 
plane y=0.1 is displayed in page 20. 

A second Taylor-Galerkin flow solution was computed by using the 
same parameters as in the first mesh. The control file is included in 
page 21. The flow solution obtained after 4000 timesteps starting from 
free stream values is displayed as follows: 
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page 22: 
page 23: 
page 24: 


page 25: 


pressure contours on the plane y=0.1 
mach number contours on the plane y=0.1 
pressure coefficient on the su rface of the 
wedge at y=0.05 compared with the exact 
solution 

pressure coefficient at the line z=0.0 y=0.1 


The new solution presents an overall improvement in definition of 

the flow features respect to the first solution as it can be observed in 

pages 26 and 27 which show the comparison between the pressure 

solutions obtained in the two meshes. Both solutions show a fair 

agreament with the exact solution. 


The convergence curves for the residuals of the derivative of the 

density for both solutions are displayed in page 28. We have divided 

the residuals by the second value of the residual (the first one is 

usually meaningless and it is ignored) to be able to plot both curves 
together. 


The slope of the residuals in the first solution is not steeper than the 
slope for the second solution as one would have expected. An 
explanation for this is the fact that the first mesh contains a distorted 
element that induces a disturbance on the flow solution. This can be 
clearly noticed on the contours of Mach number in page 17. This does 
not happen in the second mesh which presents a better convergence 
curve. 


Finally, page 27 contains a drawing of the exact position of point of 
intersection between the expansion fan and the leading edge shock. 
At this point the leading edge^ changes direction and the solution in 
the neighborhood of the intersection is changed. The "reflected 
shock" on the boundary is more likely to be a feature of the flow 
caused by the interaction than a problem of the boundary conditions. 
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